/*==================================================
Project:       Targeting Social Programs
Authors:       Diether W. Beuermann
               Bridget Hoffmann        
               Marco Stampini 
               David L. Vargas
               Diego Vera-Cossio
----------------------------------------------------
Creation Date:    Jul 2024
Modification Date:   
Do-file version:    01
References:          
Output:             
==================================================*/

/* Regular PMT model using all SISBEN varibles and SISBEN income
 This will serve as baseline for comparison		*/


/*==================================================
              0: Program set up
==================================================*/
*Written on STATA 17
drop _all
 
** source dir
cd "${dir5r}" // graph dir

/*==================================================
             1: load and transformations
==================================================*/

*----------- 1.0.1 Load CRRA welfare FUN
qui do "${dir6r}/ados/bstrap_pmt.ado" // load bstap FUN
qui do "${dir6r}/ados/bstrap_pmt_mh.ado" // load bstap FUN
qui do "${dir6r}/ados/bstrap_pmt_mh4.ado" // load bstap FUN
qui do "${dir6r}/ados/baseline_pmt.ado" // load baseline FUN
qui do "${dir6r}/ados/expanded_pmt.ado" // load expanded FUN

*----------  1.1. Data prep:
use "${dir3r}/01_survey/survey_targeting_r1_long.dta", clear

* new interactions 
cap drop d_loss_illness_inc
cap drop d_reco_ilness_inc
gen d_loss_illness_inc = illness_incidence * d_loss_ilness
gen d_reco_ilness_inc = lag_illness_incidence * d_reco_ilness
lab var d_loss_illness_inc	"Loss: Illness Incidence"
lab var d_reco_ilness_inc	"Recovery: Illness Incidence"
cap drop d_loss_cut_remittance_inc
cap drop d_reco_cut_remittance_inc
gen d_loss_cut_remittance_inc = cut_remittance_incidence * d_loss_cut_remittance
gen d_reco_cut_remittance_inc = lag_cut_remittance_incidence * d_reco_cut_remittance
lab var d_loss_cut_remittance_inc "Loss: Remittances Cut Incidence"
lab var d_reco_cut_remittance_inc "Recovery: Remittances Cut Incidence"

* Globals with shocks
glo loss d_loss_ilness d_loss_cut_remittance death divorce bankrupt nat_shock
glo loss_rec      d_loss_ilness d_reco_ilness d_loss_cut_remittance d_reco_cut_remittance death divorce bankrupt nat_shock 
glo loss_inc      illness_incidence cut_remittance_incidence death divorce bankrupt nat_shock 
glo loss_rec_inc  d_loss_illness_inc d_reco_ilness_inc d_loss_cut_remittance_inc d_reco_cut_remittance_inc death divorce bankrupt nat_shock 

*Globals with updted survey variables
glo hh_char_srv urban age_feb20_srv prop_kids_srv n_members_srv elementary_srv highschool_srv tertiary_edu_srv children_srv living_couple
glo assets_srv own_washing_machine_srv own_tractor_srv own_moto_srv own_car_srv own_PC_srv own_fridge_srv
glo hh_srv  wall_finished_srv floor_finished_srv cookpower_connected_srv wc_flush_srv kitchen_srv electricity_srv running_water_srv trash_service_srv 
glo labour_srv formal_work_srv informal_work_srv
// Including labour status from SISBEN
glo pmt_1_srv  $hh_char_srv $assets_srv $hh_srv $labour_srv
** make sure all labels are OK 
lab var illness_incidence 			"Illness Incidence"
lab var cut_remittance_incidence 	"Remittances Cut Incidence"
lab var death						"Death"
lab var divorce					    "Separation of Spouses"
lab var bankrupt					"Bankruptcy or Business Closure"
lab var nat_shock					"Natural Disaster or Fire"
lab var job_loss                   "Loss: Job lost at any point of the year"
lab var job_loss_formal            "Loss: Job lost formal at any point of the year"
lab var job_loss_informal          "Loss: Job lost informal at any point of the year"


/*==================================================
            2: Estiamtion
==================================================*/

*----------  2.1 Estimate PMT correcting for moral hazard for survey data:

// -------- Bootstrap iterations
local nreps=1000

//define poverty line in terms of Colombian Pesitos
local base=1.90*1515.1
// 2019 PPP factor https://data.worldbank.org/indicator/PA.NUS.PRVT.PP?locations=CO
 
//============== Loop around moral hazard levels

forvalues k = 1(1)5 {

   noi di "MH - inc q: `k'"

   bstrap_pmt_mh4,  incvar(d_income_nt) variables(d_gain_work_any_srv d_loss_work_any_srv) reps(`nreps')  mhazard(0.087) gainvars(d_gain_work_any_srv) lossvars(d_loss_work_any_srv) base(`base') pmtvars(l_pp_inc_srv $pmt_1_srv) predvars($pmt_1_srv) mhincq(`k')
   mat amat = (`k' \ `k' \ `k') , ( 0.087 \ 0.087 \ 0.087) , r(estt)
   if (`k' == 1)     mat estmat01 = amat
   else              mat estmat01 = estmat01 \ amat


   noi di "MH 0.15 - inc q: `k'"

   bstrap_pmt_mh4,  incvar(d_income_nt) variables(d_gain_work_any_srv d_loss_work_any_srv) reps(`nreps')  mhazard(0.15) gainvars(d_gain_work_any_srv) lossvars(d_loss_work_any_srv) base(`base') pmtvars(l_pp_inc_srv $pmt_1_srv) predvars($pmt_1_srv) mhincq(`k')
   mat amat = (`k' \ `k' \ `k'), ( 0.15 \ 0.15 \ 0.15) ,  r(estt)
   if (`k' == 1)     mat estmat02 = amat
   else              mat estmat02 = estmat02 \ amat

   noi di "MH 0.2 - inc q: `k'"

   bstrap_pmt_mh4,  incvar(d_income_nt) variables(d_gain_work_any_srv d_loss_work_any_srv) reps(`nreps')  mhazard(0.2) gainvars(d_gain_work_any_srv) lossvars(d_loss_work_any_srv) base(`base') pmtvars(l_pp_inc_srv $pmt_1_srv) predvars($pmt_1_srv) mhincq(`k')
   mat amat = (`k' \ `k' \ `k'), ( 0.2 \ 0.2 \ 0.2), r(estt)
   if (`k' == 1)     mat estmat03 = amat
   else              mat estmat03 = estmat03 \ amat
}



/*==================================================
            3: Consolidate data
==================================================*/

** send to stata
cap frame create points
frame change points 

*** consolidate matrices 

mat estmat0 = estmat01 \ estmat02 \ estmat03

** import estiamtes matrix as a dataframe
drop _all 
svmat estmat0, names(col)
gen type_m = "HM_rob"
tempfile dfdyn
save `dfdyn', replace

cap drop base _merge 
rename c1 inc_q
rename c2 mh_level

foreach v in covered inclusion_err inclusion_ok exclusion_err exclusion_ok zero epov {
    replace `v' = `v' * 100
	replace `v'_uci=`v'_uci*100
	replace `v'_lci=`v'_lci*100
}
 // PPP values 

 gen hh_benefits_lcu = hh_benefits
 replace hh_benefits = hh_benefits / 1515.1 
 replace hh_benefits_uci= hh_benefits_uci/1515.1
 replace hh_benefits_lci= hh_benefits_lci/1515.1
 // 2019 PPP factor https://data.worldbank.org/indicator/PA.NUS.PRVT.PP?locations=CO

// add previous results
append using "${dir3r}/03_auxiliar/Boostrap_estimates.dta"

foreach v in covered hh_benefits inclusion_err inclusion_ok exclusion_err exclusion_ok crra crra_l crra_u zero epov covered_uci covered_lci inclusion_err_uci inclusion_err_lci inclusion_ok_uci inclusion_ok_lci exclusion_err_uci exclusion_err_lci exclusion_ok_uci exclusion_ok_lci crra_uci crra_lci crra_l_uci crra_l_lci crra_u_uci crra_u_lci hh_benefits_uci hh_benefits_lci budget_nat_uci budget_nat_lci zero_uci zero_lci epov_uci epov_lci {
    sum `v' if year == 2019 & model_n == 1
    replace `v' = r(mean) if year == 2019 & model_n != 1 
}

// store data for later use 
save "${dir3r}/03_auxiliar/Boostrap_estimates_F5_mh.dta" , replace

/*==================================================
            4: Figure plot
==================================================*/


*--------------- Summary plot 
preserve

keep if year == 2020 
sum crra if model_n == 1 // CRRA PMT
sca crra_pmt = r(mean)

sum crra if model_n == 6 // CRRA main 
sca crra_main = r(mean)

replace inc_q = inc_q + 0.1 if mh_level == float(0.15)
replace inc_q = inc_q + 0.2 if mh_level == float(0.2)

// Exclusion
keep if year == 2020 
sum exclusion_err if model_n == 1 // CRRA PMT
sca exclusion_err_pmt = r(mean)

sum exclusion_err if model_n == 6 // CRRA main 
sca exclusion_err_main = r(mean)

#delimit ;
tw  (scatter exclusion_err inc_q if mh_level == float(0.087), mcolor(red%80))
    (rcap exclusion_err_uci exclusion_err_lci inc_q if mh_level == float(0.087), lcolor(red%30) ) 
    (scatter exclusion_err inc_q if mh_level == float(0.15), mcolor(blue%80))
    (rcap exclusion_err_uci exclusion_err_lci inc_q if mh_level == float(0.15), lcolor(blue%30) ) 
    (scatter exclusion_err inc_q if mh_level == float(0.2), mcolor(gray%80))
    (rcap exclusion_err_uci exclusion_err_lci inc_q if mh_level == float(0.2), lcolor(gray%30) ) 
    ,
    yline(`=exclusion_err_pmt', lpattern(-))  
    yline(`=exclusion_err_main', lpattern(...-) lcolor(black) )
    xtitle("Income quintile")
    ytitle("Exclusion Error")
    xlabel(1(1)5.5)
    legend(off)
name(A, replace)
;
#delimit cr

// inclusion 
keep if year == 2020 
sum inclusion_err if model_n == 1 // CRRA PMT
sca inclusion_err_pmt = r(mean)

sum inclusion_err if model_n == 6 // CRRA main 
sca inclusion_err_main = r(mean)

#delimit ;
tw  (scatter inclusion_err inc_q if mh_level == float(0.087), mcolor(red%80))
    (rcap inclusion_err_uci inclusion_err_lci inc_q if mh_level == float(0.087), lcolor(red%30) ) 
    (scatter inclusion_err inc_q if mh_level == float(0.15), mcolor(blue%80))
    (rcap inclusion_err_uci inclusion_err_lci inc_q if mh_level == float(0.15), lcolor(blue%30) ) 
    (scatter inclusion_err inc_q if mh_level == float(0.2), mcolor(gray%80))
    (rcap inclusion_err_uci inclusion_err_lci inc_q if mh_level == float(0.2), lcolor(gray%30) ) 
    ,
    yline(`=inclusion_err_pmt', lpattern(-))  
    yline(`=inclusion_err_main', lpattern(..-) lcolor(black) )
    xtitle("Income quintile")
    ytitle("Inclusion Error")
    xlabel(1(1)5.5)
    legend(off)
    name(B, replace)
;
#delimit cr

// Transfersize
keep if year == 2020 
sum hh_benefits if model_n == 1 // CRRA PMT
sca hh_benefits_pmt = r(mean)

sum hh_benefits if model_n == 6 // CRRA main 
sca hh_benefits_main = r(mean)

#delimit ; 
tw  (scatter hh_benefits inc_q if mh_level == float(0.087), mcolor(red%80))
    (rcap hh_benefits_uci hh_benefits_lci inc_q if mh_level == float(0.087), lcolor(red%30) ) 
    (scatter hh_benefits inc_q if mh_level == float(0.15), mcolor(blue%80))
    (rcap hh_benefits_uci hh_benefits_lci inc_q if mh_level == float(0.15), lcolor(blue%30) ) 
    (scatter hh_benefits inc_q if mh_level == float(0.2), mcolor(gray%80))
    (rcap hh_benefits_uci hh_benefits_lci inc_q if mh_level == float(0.2), lcolor(gray%30) ) 
    ,
    yline(`=hh_benefits_pmt', lpattern(-))  
    yline(`=hh_benefits_main', lpattern(...-) lcolor(black) )
    xtitle("Income quintile")
    xlabel(1(1)5.5)
    ytitle("Per-household Transfer ($ USD PPP)")
    legend(off)
name(C, replace)
;
#delimit cr

// rho = 3 
#delimit ;
tw  (scatter crra inc_q if mh_level == float(0.087) , mcolor(red%80))
    (rcap crra_uci crra_lci inc_q if mh_level == float(0.087), lcolor(red%30) ) 
    (scatter crra inc_q if mh_level == float(0.15), mcolor(blue%80))
    (rcap crra_uci crra_lci inc_q if mh_level == float(0.15), lcolor(blue%30) ) 
    (scatter crra inc_q if mh_level == float(0.2), mcolor(gray%80))
    (rcap crra_uci crra_lci inc_q if mh_level == float(0.2), lcolor(gray%30) ) 
    ,
    yline(`=crra_pmt', lpattern(-))  
    yline(`=crra_main', lpattern(..-) lcolor(black) )
    xtitle("Income quintile")
    ytitle("Social Welfare (CRRA: {&rho} = 3)")
    legend(order(1 "MH 8.7%" 3 "MH 15%" 5 "MH 20%") row(1))
    legend(off)
    xlabel(1(1)5.5)
    name(D, replace)
;
#delimit cr 

restore


*net install grc1leg,from( http://www.stata.com/users/vwiggins/)
grc1leg A B C D, legendfrom(D)
graph export "F5_moral_harzard_rob_inc_q.png", as(png) replace width(1200)
graph export "F5_moral_harzard_rob_inc_q.pdf", as(pdf) replace 

